clear all;clc
addpath(genpath('.'));  %添加工具箱路径
% load Stego,11-3号来测试image3
%% 超参数设置
Block_Size=128;%检测窗口大小s
k=64/256;%前面测量的步长都是尽可能的小的，原文是Block_sie的1/4，我这里是15/64
step=fix(k*Block_Size);%当降低为32，步长只有7，还是向下取整，初值7.5
save step step
%% 读入待检测的屏摄水印图像，晚上测image5，侧拍  
% Stego=rgb2ycbcr(imread('./Ext_exp/Fang/image1/A30_6.jpg'));
Stego=rgb2ycbcr(imread('./distance_exp/image5/Cloc_1_dc.jpg'));
VAR='5Cloc_1_dc';
%% 多尺度响应值计算   
Y_Stego_A=Stego(:,:,1);  
%% 上采样
% Y_Stego_A_1=imresize(Y_Stego_A,1.1);
% Y_Stego_A_2=imresize(Y_Stego_A,1.2);
% Y_Stego_A_3=imresize(Y_Stego_A,1.3);
% Y_Stego_A_4=imresize(Y_Stego_A,1.4);
%% 下采样
Y_Stego_B=imresize(Y_Stego_A,0.9);
Y_Stego_C=imresize(Y_Stego_A,0.8);
Y_Stego_D=imresize(Y_Stego_A,0.7);
Y_Stego_E=imresize(Y_Stego_A,0.6);
Y_Stego_F=imresize(Y_Stego_A,0.5);
Y_Stego_G=imresize(Y_Stego_A,0.4);
Y_Stego_H=imresize(Y_Stego_A,0.3);
Y_Stego_I=imresize(Y_Stego_A,0.2);
% Y_Stego_J=imresize(Y_Stego_A,0.1);
% Y_Stego_K=imresize(Y_Stego_A,0.05);
Y_Stego_Set={Y_Stego_A,Y_Stego_B,Y_Stego_C,Y_Stego_D,Y_Stego_E,Y_Stego_F,Y_Stego_G,Y_Stego_H,Y_Stego_I};%,Y_Stego_J,Y_Stego_K};
% 针对上采样
% Y_Stego_Set={Y_Stego_A_4,Y_Stego_A_3,Y_Stego_A_2,Y_Stego_A_1,Y_Stego_A,Y_Stego_B,Y_Stego_C,Y_Stego_D,Y_Stego_E,Y_Stego_F,Y_Stego_G,Y_Stego_H,Y_Stego_I};%,Y_Stego_J,Y_Stego_K};
R_Scalar={};
Scale=[1,0.9,0.8,0.7,0.6,0.5,0.4,0.3,0.2];%,0.1,0.05];
% Scale=[1.4,1.3,1.2,1.1,1,0.9,0.8,0.7,0.6,0.5,0.4,0.3,0.2];%,0.1,0.05];
R_Intensy=zeros(size(Scale));
count=1;
%% 依次计算每个尺度下的响应值
for Y_Stego=Y_Stego_Set
    Y_Stego=Y_Stego{1};  
%     Y_Stego=Y_Stego_G;%测试一下,并且分析一下
    %加上全图直方图纹理增强，不好用，一般是先图像增强再维纳滤波
%     Y_Stego=adapthisteq(Y_Stego,'clipLimit',1,'NumTiles',[32 32]);%自适应直方图均衡纹理增强
    [w,h]=size(Y_Stego);
    % 全局维纳滤波默认3*3窗口，可用
    Y_Stego=Y_Stego-wiener2(Y_Stego);
%     J=histeq(Y_Stego,64);
%     row_var=var(double(Y_Stego));col_var=var(double(Y_Stego'));
%     row_mean=sum(Y_Stego,1);col_num=sum(Y_Stego,2);
%     [~,LOCS_row]=findpeaks(row_num);[~,LOCS_col]=findpeaks(col_num);
    Y_Stego=My_padding(Y_Stego,Block_Size);
    for i=1:step:(w)
        for j=1:step:(h)
            Detect_Block=Y_Stego((i:i-1+Block_Size),(j:j-1+Block_Size));
            %局部直方图增强+维纳滤波
%             Detect_Block=adapthisteq(Detect_Block,'clipLimit',1,'NumTiles',[32 32]);%自适应直方图均衡纹理增强
%             Detect_Block=Detect_Block-wiener2(Detect_Block);
            [DIFF((i-1)/step+1,(j-1)/step+1),R((i-1)/step+1,(j-1)/step+1)]=My_Bline_Detector(Detect_Block);%R中存储了每个窗口的同步响应值
        end
    end
    R_Scalar{count}=R;
    %% 标注一下最大异常对应的块
%     figure,
%     plot_area(R,Y_Stego,Block_Size);
    R_Intensy(count)=std2(R);%*Scale(count);%
    R=[];
    count=count+1;
end

%% 依据强度进行尺度选择
% [~,id]=max(R_Intensy);
% disp(['最终尺度选择:',num2str(Scale(id))])
%% 依据峰值进行尺度选择
[~,id]=findpeaks(R_Intensy);%,'minpeakheight',mean(C_row));
if isempty(id)==1
    id=1;
end
disp(['最终尺度选择:',num2str(Scale(id(1)))])
id=id(1);
%% 3*3窗口维纳滤波获取噪声分量
% J = imnoise(Y_Cover,'gaussian',0,0.005);
save id id
save Scale Scale
save R_Scalar R_Scalar

save([VAR '.mat'], 'DIFF');